## This program generates Figure 1 from "Financial Constraints, 
##  Sectoral Heterogeneity, and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(sandwich)
library(tidyverse)
library(glue)
library(dynlm)
library(seasonal)

## Set maximum length of response horizon (default:16)
numres <- 16 

## Set number of standard deviations to show in errors (default: 1.65)
se_scale <- 1.65

## Set number of lags of the dependent variable and shock (default: 8)
lags <- 8

## Specify regression formula
f <- formula( glue(paste("L(y,-i) ~ trend(y) + season(y) + mp_shock", 
                        "+ L(mp_shock,1:lags)  ", 
                        " +L(y,1:lags)" )))

## List of different MP shocks to use
shocknames <- c("rr","gk","mar","brw")

## Full names of each shock to use in plots
longshocknames <- c("Romer and Romer (2004)", "Gertler and Karadi (2015)",
                    "Miranda-Agrippino and Ricco (2021)", "Bu, Rogers, and Wu (2021)")

## Beginning dates will vary across series since they cover different periods, but
##  all the series will end in 2008Q4.
fin   <- c(2008,4) ## Default: c(2008,4)
                

## Load data 
data <- ts(read.csv("Aggregate_data.csv"),start=c(1966,1),end=c(2021,4),frequency=4)

# Shock coefficient estimates and standard errors
gamma    <- matrix(0,nrow=(numres+1),ncol=length(shocknames))
gamma_se <- matrix(0,nrow=(numres+1),ncol=length(shocknames))




##### Step 1: Run regressions #####



## Extract the log of the real aggregate NPPE series
y <- log(data[,"rnppe_all"])


## Outer loop over each type of shock
for(j in 1:length(shocknames)){
  
  ## Specify the monetary policy shock
  eval(parse(text=paste0("mp_shock <- data[,'",shocknames[j],"_shock']")))
  
  ## Loop over response horizon
  for(i in 0:numres){
    
    ## Each shock starts at different times
    if(j==1){
      
      ## rr shocks start in 1970Q1
      begin <- c(1970,1)
      
    } else if(j==2){
      
      ## gk shocks start in 1975Q1
      begin <- c(1975,1)
      
    } else if(j==3){
      
      ## mar shocks start in 1991Q1
      begin <- c(1991,1)
      
    } else {
      
      ## brw shocks start in 1995Q1
      begin <- c(1995,1)
      
    }
    
    ## Need to calculate date of last shock used so all data end in 2008. For
    ##  example, if the data end in 2008, then the latest shock used when estimating
    ##  the response 13 quarters after the shock will need to be in the third
    ##  quarter of 2005, since (2005Q3 + 13 quarters = 2008Q4).
    lastshock_y <- fin[1]-floor((i+3)/4)
    lastshock_q <- 4-i%%4
    
    ## Estimate regression
    eq <- dynlm(formula=f,start=begin,end=c(lastshock_y,lastshock_q),data=data)
    newey <- NeweyWest(eq,prewhite=FALSE)
    
    ## Save coefficient estimates
    gamma[(1+i),j] <- coef(summary(eq))["mp_shock","Estimate"]
    gamma_se[(1+i),j] <- sqrt(newey["mp_shock","mp_shock"])
    
  }
  
}



##### Step 2: Plot IRFs #####

## Set scale of each type of shock 
## rr: 1pp
## gk: 1pp
## mar: 0.054 (two standard deviations of non-zero values)
## brw: 0.12 (two standard deviations of non-zero values)

shock_scale <- c(1,1,0.054,0.12)

qtrs <- c(0:numres) ## Generates x-axis sequence for plots

## Outer loop over shock types
for(j in 1:length(shocknames)){
  
  ## Create new window
  windows(width=5,height=4)
  par(oma=c(2,2,1,0))
  
  ## Specify which scale to show for each shock
  scale <- shock_scale[j]
  
  irf_agg1 <- scale*(gamma[,j])
  irf_agg2 <- scale*(gamma[,j] + se_scale*gamma_se[,j])
  irf_agg3 <- scale*(gamma[,j] - se_scale*gamma_se[,j])
  
  yrange <- c(-0.04,0.05)
  
  # Plot response of total
  par(mar=c(2,2,2,1))
  plot(x=qtrs,y=irf_agg1,type="b",col="blue",lwd=3,ylim=yrange,pch=1,
       axes=F,frame=T,ann=FALSE,cex.main=1,cex=1.5)
  axis(1,cex.axis=1.75,at=pretty(qtrs),lab=pretty(qtrs))
  axis(2,cex.axis=1.75,at=pretty(yrange),lab=paste0(pretty(yrange)*100, "%"),las=TRUE)
  {mtext("Quarters after shock",side=1,line=2.5,cex=1.5)}
  abline(a=0,b=0,col="black",lwd=1)
  
  lines(x=qtrs,y=irf_agg2,lty=2,col="red",lwd=2)
  lines(x=qtrs,y=irf_agg3,lty=2,col="red",lwd=2)
  
  title(main=longshocknames[j], col.main="blue", font.main=2,cex.main=1.25)
  
  
}
